Industrial 

Mathematic 

Institute 


1997:16 

A  characteristic  domain 
decomposition  technique  for  two- 
phase  flows  with  interfaces 


H.  Wang  and  B.G.  Ersland 


1 


A  Characteristic  Domain 
Decomposition  Algorithm  for 
Two-Phase  Flows  with  Interfaces 

HONG  WANG  and  BRIT  GUNN  ERSLAND 


1.1  INTRODUCTION 

Tlie  mathematical  model  that  describes  the  process  of  an  immiscible  displacement  of 
oil  by  water  in  reservoir  production  or  other  two-phase  fluid  flows  in  porous  media  leads 
to  a  strongly  coupled  system  of  a  degenerated  nonlinear  advection-diffusion  equation 
for  saturation  and  an  elliptic  equation  for  pressure  and  velocity.  The  hyperbolic 
nature,  strong  coupling,  and  nonlinearity  of  the  system  and  the  degeneracy  of  the 
diffusion  makes  numerical  simulation  a  challenging  task.  Many  numerical  methods 
suffer  from  serious  non-physical  oscillations,  excessive  numerical  dispersion,  and/or 
a  combination  of  both  [CJ86,  Ewi84].  Previously,  Espedal,  Ewing,  and  coworkers 
developed  a  characteristic,  operator-splitting  technique  in  effectively  solving  two-phase 
fluid  flow  problems  [DEES90,  EE87].  In  practice,  a  reservoir  often  consists  of  different 
subdomains  with  different  porosities  and  permeabilities.  In  the  case  of  single-phase 
fluid  flows  the  concentration  and  total  flux  are  continuous  across  the  interfaces  between 
different  subdomains  since  the  diffusion  never  vanishes.  Our  earlier  work  addressed 
numerical  simulation  to  linear  transport  equations  arising  in  single-phase  flows  with 
interfaces  [?].  However,  in  the  case  of  two-phase  fluid  flows  the  saturation  equation 
itself  is  nonlinear  and  different  subdomains  have  different  capillary  pressure  curves. 
The  continuity  of  capillary  pressure  across  interfaces  implies  a  jump  discontinuity  of 
the  water  saturation  at  the  same  locations.  The  jump  discontinuity  of  the  saturation 
at  the  interfaces  might  incur  some  oscillations  around  the  interfaces,  which  can  be 
propagated  into  the  domain  and  destroy  the  overall  accuracy.  Hence,  great  care  has 
to  be  taken  in  the  development  of  an  effective  solution  procedure  for  the  simulation 
of  two-phase  fluid  flows  in  porous  media  with  interfaces. 

This  paper  describes  a  characteristic-based,  non-overlapping  domain  decomposition 
algorithm  for  solving  the  saturation  equation  in  two-phase  fluid  flows  with  interfaces. 
First,  with  the  known  saturation  at  the  previous  time  step  one  obtains  an  approximate 
Dirichlet  boundary  condition  at  the  outflow  domain  interface  by  integrating  the 
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saturation  equation  (ignoring  the  capillary  pressure  term)  along  characteristics.  With 
the  approximate  outflow  Dirichlet  boundary  condition  at  the  domain  interface  and  the 
given  boundary  condition  at  the  physical  inflow  boundary  one  can  close  the  system 
on  the  current  subdomain  and  applies  the  characteristic  operator-splitting  procedure 
[DEES90.  EE87]  to  solve  the  full  saturation  equation  (including  the  capillary  pressure 
effect).  Second,  one  uses  the  continuity  of  capillary  pressure  across  the  domain  interface 
to  pass  the  value  of  saturation  as  an  approximate  inflow  Dirichlet  boundary  condition 
to  the  next  subdomain,  one  then  applies  the  same  characteristic  operator-splitting 
procedure  to  solve  the  saturation  equation  on  the  current  subdomain.  Third,  according 
to  the  overall  loss  or  gain  of  mass  one  adjusts  the  approximate  outflow  Dirichlet 
boundary  condition  at  the  domain  interface  to  iterate  between  different  subdomains 
until  the  algorithm  converges.  Finally,  a  mixed  method  is  adopted  to  solve  the  pressure 
equation  due  to  its  accurate  approximation  to  the  velocity  held  and  its  local  mass 
conservation  property. 

The  rest  of  the  paper  is  organized  as  follows:  In  Sections  2  and  3  we  formulate  the 
problem  and  discuss  related  solution  techniques.  In  Section  4  we  present  a  domain 
decomposition  algorithm  for  the  two-pliase  fluid  how  problems  with  interfaces.  In 
Section  5  we  present  some  numerical  results  to  show  the  promise  of  the  method. 


1.2  PROBLEM  FORMULATION 

A  suitable  mathematical  model  for  the  total  Darcy  velocity  m,  the  total  pressure  p,  and 
the  water  saturation  S  G  [0, 1]  in  an  incompressible  displacement  of  oil  by  water  in  a 
porous  medium  can  be  described  by  the  following  set  of  partial  differential  equations 

[CJ86]: 


(x,t)  G  SI  x  [0,T], 

(x.l)  f  !!x  [0.  Y'\  (1) 

(x,t)  G  Oil  x  [0,T], 


(x ,t)  g!!x  [0,T], 

(x,  t)  G  30  x  [0.  T], 
x  G  ft, 

where  ft  is  the  physical  domain,  K_(x)  is  the  absolute  permeability  tensor  of  the 
medium,  A*,  i  =  o,w,  denotes  the  water  and  oil  mobilities  respectively,  qi(x,t.)  and 
qs(x,t)  are  source  terms,  f/zlift)  and  q4(x,t)  are  the  prescribed  boundary  conditions, 
n(x)  is  the  unit  outward  normal  vector,  <<  1  is  a  scaling  factor  to  the  diffusion  term, 
pc  is  the  capillary  pressure,  and  f(S)  and  D_(S,x)  are  the  fractional  how  function  and 
diffusion  term  given  by 


V  •  -u(x,  i)  =  qi(x,  t), 

u(x,  t)  =  -K(x)(X0(S)  +  \w(S))Vp(x,  t), 
u(x,  t)  ■  n(x)  =  q2 §1'  t) i 

and 

r\  q 

<Kx)—  +  V  ■  (f(S)u- eD(S,x)VS)  =  /  ). 

(f(S)u  —  eD_(S,  x)V S)  ■  n(x)  =  q4($,  t), 
S(x.  0)  =  Sq{x), 


f(S) 

D(S.x) 


K(S) 

K(S)  +  A0(S) 

XW(S)X0(S)  dPc 

-(-}  \0(S)  +  \W(S)  dS' 


(3) 
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Note  that  the  two  equations  in  (1)  form  a  second-order  elliptic  equation  for  the 
pressure  p(x,  t)  and  are  coupled  to  the  saturation  equation  (2)  through  the  saturation 
5  in  the  coefficients.  On  the  other  hand,  saturation  equation  (2)  is  a  nonlinear 
advection-diffusion  equation  and  is  coupled  to  the  pressure  equation  (1)  through  the 
Darcy  velocity  u.  Furthermore,  in  the  mathematical  model  the  diffusion  term  D(  S.x) 
vanishes  at  5  =  0  and  5=1,  which  is  an  idealized  case  since  physically  D( S.x) 
vanishes  for  S  G  [0,  5rr]  or  S  G  [1  —  S" .  1]  with  S "  being  the  irreducible  saturation 
value.  The  fractional  flow  function  f(S)  defined  in  (3)  is  typically  an  5-shaped  curve 
of  saturation  5  and  degenerates  at  5  =  0  (with  the  same  understanding).  Because  the 
saturation  profile  is  usually  a  decreasing  function  in  space,  as  time  evolves  /(5)  tends 
to  force  a  shock  discontinuity  to  develop  in  5  while  the  diffusion  term  D_(S.x)  tends 
to  prevent  the  shock  from  forming.  The  dynamic  process  could  be  fairly  complicated 
because  the  diffusion  degenerates  in  front  of  the  steep  saturation  front.  It  depends  on 
the  interaction  between  the  advection  and  diffusion  terms,  in  particular,  on  the  rates 
at  which  D_(S.x)  and  f(S)  tend  to  zero  as  5  tends  to  zero. 

When  the  physical  domain  O  is  composed  of  different  media,  the  different  porosities 
and  permeabilities  result  in  different  capillary  pressure  curves  on  each  subdomain 
(Figure  ??).  Across  an  interface  T  the  phase  pressures  are  continuous  and  mass  is 
conserved,  leading  to  the  following  interface  conditions 

Pc(S)  |r_  =Pc(S)  |r+, 

u-n|r_  =u-n|r+,  (4) 

(f(S)u  -  eD(S ,  s)V5)  ■  njr_  =  (f(S)u  -  eD{S ,  x)V  S)  ■  n|r+  - 

The  continuity  of  capillary  pressure  pc  across  an  interface  F  implies  the  discontinuity 
of  the  saturation  across  the  interface  (Figure  ??).  One  has  to  resolve  the  discontinuity 
carefully  so  that  no  spurious  effects  will  be  propagated  into  the  domain. 


1.3  OPERATOR  SPLITTING  TECHNIQUES 

Extensive  research  has  been  carried  out  for  the  numerical  simulation  of  system  (  I  I  (2) 
without  interfaces.  Various  techniques  have  been  developed  to  decouple  and  linearize 
the  system,  including  a  fully  coupled  and  fully  implicit  linearization  strategy,  an 
IMPES  (IMplicitly  advances  the  Pressure  and  Explicitly  updates  the  Saturation  in 
time)  strategy,  and  a  sequential  time  stepping  strategy  [Ewi84].  Different  numerical 
methods,  including  the  standard  Galerkin  finite  element  method,  the  cell-centered 
finite  difference  method,  the  finite  volume  method,  and  the  mixed  finite  element 
method,  have  been  used  to  solve  the  pressure  equation  [CJ86,  DEES90.  DEW83,  EE87, 
Ewi84].  We  used  the  mixed  method  to  solve  the  pressure  equation  due  to  its  accurate 
approximation  to  the  velocity  held  and  its  local  mass  conservation  property.  Because 
the  normal  component  of  the  velocity  held  is  continuous,  the  discrete  algebraic  system 
for  the  pressure  equation  is  in  fact  the  same  as  that  with  no  interfaces.  Hence,  one  can 
solve  the  global  system  as  usual.  Alternatively,  one  can  use  a  domain  decomposition 
procedure  to  solve  the  pressure  equation  on  each  subdomain  iteratively.  We  refer  the 
interested  readers  to  [BW86,  SBG96]  and  the  references  therein  for  details. 

For  simplicity  of  exposition  we  consider  a  one-dimensional  analogue  of  equation  (2). 
Notice  that  equation  (2)  is  almost  hyperbolic  due  to  the  small  parameter  e  «  1.  An 


4 


H.  WANG  AND  B.G.  ERSLAND 


effective  solution  procedure  for  solving  the  dominating  adveetive  part  of  equation  (2) 

*(x)^  +  -k{f{S)u)  =  0  (5) 

is  to  discretize  equation  (5)  along  the  characteristics,  which  allows  large  time  steps 
to  be  used  in  the  numerical  simulation.  Because  equation  (5)  may  have  more  than 
one  solution  due  to  the  shape  of  the  fractional  flow  function  f(S),  one  cannot  directly 
apply  the  modified  method  of  characteristics  [DR82]  to  equation  (5).  We  follow  the 
work  of  Espedal,  Ewing,  and  coworkers  [DEES90,  EE87]  and  split  the  fractional  flow 
function  f(S)  into  two  parts  by 


with 


f(S)  =  f(S)  +  b(S)S , 
/( Sbl ) 


m  = 


Sbl 

f(S). 


-S,  ifO  <S<Sbl, 


if  Sbl  <  S  <1. 
Here  the  Buckley-Leverett  shock  saturation  Sbl  is  defined  by 

/( Sbl ) 


/  (Sbl)  — 


'  BL 


(6) 

(7) 

(8) 


Because  / ( S)u  gives  the  unique  physical  velocity  for  an  established  shock,  we  use 
this  operator  splitting  and  rewrite  equation  (2)  along  the  characteristics  as 


flC'ii+l  ocn+l  ,  ocn+1 

=  <K'X)-^-  +  f  (Sn+1)u 


Or 


"+1  _ _  q_ 


dx 


and 


dSn+l  g 

<P(x)—K -  +  7- 

OT  OX 


/  _  _  ocn+l. 

)Sn+1u  -  eD(Sn+%,x)  j  =  q3 ps,  tn+1 ). 


(9) 


(10) 


From  the  definition  of  /  it  follows  that  the  characteristic  direction  is  uniquely 
determined  by  equation  (9)  since  the  shape  of  /  allows  only  a  rarefaction  wave  and 
a  contact  discontinuity  for  a  non-increasing  saturation  profile.  Thus,  the  hyperbolic 
equation  (9)  is  discretized  by  integrating  backwards  along  the  characteristics 


s*  =  x  -  f  [S'"  i  At. 


(ID 


where  Sn*  =  S(x*,tn)  and  At  =  tn+1  —  tn  is  the  time  step. 

Note  that  the  characteristics  determined  by  equation  (9)  are  all  straight  lines  in  the 
(a;,  t)  plane.  If  equation  (9)  is  solved  exactly,  the  only  change  in  the  solution  along  the 
characteristics  is  due  to  diffusion  (and  possibly  the  source  term  which  vanishes  except 
at  wells).  Thus,  we  solve  equation  (10)  by  the  modified  method  of  characteristics 
[DEES90,  DR82.  EE87] 


/ 

JQ 


£fn-\- 1  gn 

At. 


-wdfl  + 


Jn 


Vw(x)  €  Hq  (fi). 


J  OX 


(12) 
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Hero  a  characteristic  tracking  is  used  for  the  advection  term,  and  the  quadratic  Petrov- 
Galerkin  method  is  used  for  the  diffusion  term  and  the  residual  advection  term  where 
the  trial  functions  S  are  chosen  to  be  hat  functions  and  the  test  functions  w(x)  are 
constructed  by  adding  an  quadratic  perturbation  to  the  hat  functions  [DEES90,  EE87]. 


1.4  A  CHARACTERISTIC  DOMAIN  DECOMPOSITION 

ALGORITHM  FOR  SYSTEM  (1)  (2)  WITH  INTERFACES 

We  now  describe  a  characteristic  domain  decomposition  algorithm  for  solving  the 
system  (l)-(2)  with  interfaces.  We  adopt  a  sequential  solution  strategy  to  decouple 
and  linearize  the  system  [DEES90,  DEW83].  For  the  domain  decomposition  techniques 
for  pressure  equations  with  interfaces  we  refer  the  interested  readers  to  [BW86,  SBG96] 
for  details.  We  present  the  algorithm  for  a  one-dimensional  problem  on  Q.  =  (a,  b )  with 
one  interface  at  a  <  d.  <  b.  Let  N  be  a  positive  integer,  At  =  T/N ,  and  tn  =  nAt. 

Initialization 

Substitute  the  initial  condition  5(a:,0)  for  S  in  (1)  and  solve  equations 
(1)  at  t°  by  the  mixed  method  to  obtain  the  Darcy  velocity  «n(.r),  where 
un(x)  =  u(x,  tn). 

for  n  =  0, 1 . N  —  1  do 

for  l  =  0, 1 . Im  ~  1  do 

LI.  For  l  =  0,  in  equation  (2)  approximate  un+1(x)  by  Uq+1(#|  =  un(x)  or 
2m" (as)  —  «n-1(x).  For  l  >  1,  substitute  for  S  in  (1)  and  solve 

equations  (1)  at  tn+1  by  the  mixed  method  to  obtain  the  Darcy  velocity 

«?+1- 

L2.  For  l  =  0,,  assign  £« ^(d-)  =  Sn(d*),  where  <S,”j|'1(d_)  = 

lim  Si'k(x,tn+1 )  and  d*  is  defined  in  equation  (11)  with  x  being 

x—td,x<id 

replaced  by  d.  For  l  >  1,  assign  =  S^+^kM(d^). 

L3.  Use  the  interface  condition  (SVt1  (d,-))  =  p^iSVt1  (d+))  to  evaluate 
S'^1(d+),  where  -SV.V’fd.,  =  '  lim  S^(x).  ' 

x^-d,x>d  ’ 

for  k  =  0, 1 . I'm  —  1  do 

if  ERROR  >  TOLERANCE  then 

Kl.  With  the  given  inflow  boundary  condition  at  x  =  a  and 

as  the  outflow  Dirichlet  boundary  condition  at  x  =  d.  solve  equation 
(12)  on  the  subdomain  (a,d)  for  at  time  tn+1. 

K2.  with  .si;  j1  ui ..  )  as  the  inflow  Dirichlet  boundary  condition  at  x  =  d 
and  the  given  outflow  boundary  condition  at  x  =  b,  solve  equation 
(12)  on  (d.b)  for  at  tn+1  in  parallel  to  the  previous  step. 

K3.  Calculate  the  mass  error  =  AM  —  —  Sn)dfl.  where 

AM  is  the  mass  injected  at  the  inflow  boundary  and  through  the 
wells  during  the  time  period  [tn,tn+1]. 
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K4.  Update  the  Dirichlet  boundary  condition  at  the  interface  x  =  d  by 
STjHd-)  :=  S$Hd-)  +  where  K  is  a  relaxation  parameter. 

K5.  Use  the  interface  condition  (S"^1  (d-))  =  (d+) |  to 

evaluate  (d,+  ). 
else 

h  —  k,M  and  l  —  Im 

endif 

end 

k  =  Um 

end 


l  —  Im 

■un+1  :=  uV+1 

i-M 


and  Sn+1  :=  S?+\ 

Im 


end 


Note  that  the  full  equation  (12)  is  almost  symmetrized  and  almost  well  conditioned. 
Namely,  the  condition  number  is  of  order  0(DAt/(Ax)2).  Hence,  a  diagonal 
preconditioner  works  well  in  practice,  in  contrast  to  elliptic  equations  where  the 
coefficient  matrix  is  ill  conditioned  and  extensive  research  has  been  carried  out  to 
develop  an  efficient  preconditioner. 

We  now  outline  generalizations  of  the  above  algorithm  in  several  directions.  First, 
it  is  easy  to  see  that  the  above  algorithm  applies  to  problems  with  several  interfaces. 
Second,  we  note  that  the  procedure  applies  to  multidimensional  problems,  as  long  as 
the  adjustment  in  Step  K3  is  kept  local  in  space  to  avoid  introducing  any  spurious 
nonzero  saturation  to  the  location  where  the  saturation  is  zero.  Third,  Because  the 
coefficient  matrix  for  the  pressure  equation  has  a  much  bigger  condition  number  than 
that  for  the  saturation  equation,  it  is  much  more  expensive  to  solve  equations  in  (1) 
than  to  solve  equation  (2)  at  each  time  step.  Physically  the  Darcy  velocity  is  much 
smoother  and  varies  less  rapidly  than  the  saturation.  Thus,  we  can  use  larger  time 
steps  for  pressure  equations  in  (1)  and  smaller  time  steps  for  the  saturation  equation 
(2)  (see  [DEW83,  Ewi84]  for  details). 


1.5  Numerical  Experiments 


In  this  section  we  present  a  numerical  example  to  show  the  promise  of  the  algorithm. 
More  extensive  results  can  be  found  in  [Ers96].  In  the  example,  the  space  domain 
(a,  b)  =  (0,1)  with  the  interface  located  at  d  =  0.5.  The  time  interval  [0,  T]  =  [0.  0.048] , 
e  =  0.01,  A W(S)  =  S 2.  A 0(S)  =  (l-S)2.  Ax  =  1/150.  At  =  0.001.  K  =  10  on  (0,0.5) 
and  1  on  (0.5, 1).  The  initial  condition  is  an  established  shock  given  by 


Sq(x) 


0  3 

1  -  —a-.  if  0  <  x  <  0.4, 
0.4  “  “ 

0,  if  0.4  <  x  <  1. 


(13) 


In  the  numerical  experiments,  lM  =  1  and  Um  =  4.  Namely,  we  extrapolated  the 
current  velocity  field  un+1  by  its  values  at  the  previous  time  steps  and  did  not  iterate 
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on  equations  in  (1).  With  the  extrapolated  velocity  field  at  the  current  time  step,  we 
iterated  four  times  on  the  saturation  equation  (2)  at  each  time  step.  It  was  seen  in 
Figure  (??)  that  the  permeability  has  considerable  effect  on  diffusion  and  capillary 
pressure.  For  a  fixed  saturation  the  capillary  pressure  is  higher  in  a  lower  permeable 
zone  than  it  is  in  a  high  permeable  zone.  We  observe  that  the  continuity  of  capillary 
pressure  in  (4)  enforces  a  jump  up  in  the  saturation  profile  across  the  interface.  The 
numerical  results  are  free  of  oscillation  or  numerical  dispersion,  and  agree  with  the 
results  in  [CY92]. 
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